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NATIONAL  ADVISORY  COMMITTEE  FOR  AERONAUTICS 
TECHNICAL  NOTE  2539 

APPLICATION  OF  VARIATIONAL  METHODS  TO  TRANSONIC  FLOWS 
WITH  SHOCK  WAVES 

By  Chi-Teh  Wang  and  Pei-Chi  Chou 

SUMMARY 


Variational  methods  for  the  approximate  solution  of  subsonic  and 
transonic  flows  of  a  compressible  fluid  before  the  occurrence  of  shock 
waves  have  been  carried  out  in  previous  papers ♦  The  methods  fail  as 
soon’ as  the  shock  waves  occur  as  the  flow  behind  the  shock  waves  now 
becomes  rotational  and  has 'variable  entropy.  Since  most  transonic  flows 
are  accompanied  by  shock  waves ^  a  method  which  allows  for  shock  waves 
and  variable  entropy  is  necessary  for  the  study  of  such  flows.  By 
modifying  Bateman’s  variational  principle  for  irrotational  flows,  it  is 
shown  that  a  variational  principle  for  flows  with  rotation  and  variable 
entropy  can  be  obtained.  By  applying  this  variational  principle  to  the’ 
regions  of  flow  behind  shock  waves  and  Bateman’s  original  principle  to 
the  other  regions  in  the  fluid,  shock  equations  can  be  directly  obtained. 

A  procedure  for  computing  numerical  solutions  for  such  flows  is  suggested, 
and  a  numerical  example  is  carried  out.  At  high  Mach  number  above  a 
certain  limiting  value,  the  results  show  that  irrotatiorial  flow  fails. 
However,  by  inserting  shock  waves  and  allowing  a  part  of  the  flow  to  be 
rotational,  computation  indicates  that  solution  exists  again. 

INTRODUCTION  .  . 

In  previous  papers  (references  1  to  7)  the  senior  author  and  his 
associates  have  succeeded  in  applying  the  variational  methods  to  the 
study  of  subsonic  and  transonic  flows  of  a  compressible  fluid  past 
arbitrary  bodies  before  the  occurrence  of  shock  waves.  Numerical  examples 
were  carried  out  in  the  case  of  the  flow  past  a  circular  cylinder,  an  , 
elliptical  cylinder,  a  Kaplan  bump,  a  sphere,  and  an  ellipsoid.  The 
results  were  found  to  check  excellently  with  those  computed  by  other 
methods.  These  results  Indicate  that  the  variational  method  will  give 
good  approximations  to  flows  past  either  thick  or  thin  bodies  and  at 
both  low  and  high  Mach  numbers.  The  method  as  formulated,  however,  can 
only  be  .applied  to  irrotational  flows.  As  soon  as  shock  waves  occur  the 
method  fails  because  the  flow  behind  the  shocks  then  becomes  rotational. 

As  most  transonic  flows  are  accompanied  by  shock  waves,  different  methods 
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which  allow  the  rotationality  of  the  flow  as  well  as  variable  entropy 

must  be  formulated  if  any  significant  results  are  to  be  obtained  from  ^  ir 

the  studies  of  such  flows. 

A  reexamination  of  Bateman’s  variational  principles  (reference  8) 
indicates  that  with  some  modification  the  principle  in  terms  of  the 
stream  function  which  was  originally  formulated  for  irrotational 

flows,  can  also  be  applied  to  rotational  flows  with  entropy  change.  The 
resulting  variational  integral  is  the  same  as  the  one  recently  obtained 
by  Lin  and  Rubinov  (reference  9)-  In  the  study  of  transonic  flows  this 
integral  is  to  be  applied  in  the  region  of  flow  after  shock  waves  and 
the  original  Bateman  integral  in  the  other  regions.  It  can  be  shown 
that  the  shock  equations  are  directly  obtainable  from  these  principles. 

With  the  variational  principles  obtained,  a  direct  method  for  the  approxi¬ 
mate  solution  of  transonic  flow  with  shock  waves  may  again  be  formulated 
following  the  Rayleigh-Ritz  procedure.  The  actual  carrying  out  of  such 
a  method  however  was  found  too  laborious.  Instead,  Galerkin’s  method 
may  be  used  which  shortens  the  numerical  work  to  a  great  extent.  In 
this  report  approximate  solution  in  the  case  of  the  transonic  flow  past 
a  circular  cylinder  has  been  carried  out.  The  results  show  that  when 
the  Mach  number  increases  to  a  certain  limiting  value  without  allowing 
for  shock  waves  the  variational  method  does  not  have  a  solution.  This 
probably  indicates  the  breakdown  of  the  irrotational  flow.  By  allowing 
the  occurrence  of  shock  waves,  solution  again  exists. 

This  work  was  conducted  at  the  Daniel  Guggenheim  School  of  Aeronautics, 
College  of  Engineering,  New  York  University,  under  the  sponsorship  and 
with  the  financial  assistance  of  the  National  Advisory  Committee  for 
Aeronautics.  The  authors  are  indebted  to  Professor  F.  K.  Teichmann  for  ^ 

his  kind  interest. 


SYMBOLS 


a  velocity  of  sound 

undetermined  parameters 


Cp,Cv 


ci^C2,C3 

M 

n 


specific  heat  at  constant  pressure  and  volume, 
respectively 

constants 

Mach  nvmber  (q/a) 
normal  distance 
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P 

S 

u 

V 

qp 

p 

CD 

U 

Subscripts 
«  x,y,r,0 

s, 

o 


pressure  * 

magnitude  of  velocity  vector^  with  components  u 
and  V 

entropy 

velocity  component  in  x-direction 
velocity  component  in  y-direction 
coordinates  of  a  point  in  the  fluid 
isentropic  exponent  c^ 

velocity  potential  function 

stream  function 

mass  density 

rotation  of  flow 

velocity  of  undisturbed  stream 

denote  differentiation  in  corresponding  direction 

isentropic  stagnation  conditions  for  undisturbed 
stream 

conditions  of  undisturbed  stream 


GOVERNING  DIFFERENTIAL  EQUATIONS  AND  BOUNDARY  CONDITIONS 


Consider  a  steady  two-dimensional  flow  passing  a  cylinder  of 
arbitrary  shape.  Most  aerodynamic  problems  are  concerned  with  either 
flow  from  rest  or  flow  which  is  parallel  and  uniform  at  infinity.  In 
such  cases  before  the  occurrence  of  shock  waves  the  flow  is  irrotational. 
As  the  speed  of  the  flow  increases^  shock  waves  occur. ,  For  example^  the 
flow  past  an  airfoil  is  shown  in  figure  1.  The  flow  is  assumed  to  be 
subsonic  at  infinity.  Because  of  the  presence  of  the  airfoil^  if  the 
Mach  number  of  the  flow  is  sufficiently  high  there  will  be  a  region  of 
supersonic  flow  near  the  surface  of  the  airfoil^  as  shown  by  the  dotted 
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lines  in  the  figure.  Shock  waves  will  occur  at  some  points  on  the 
surface  in  these  supersonic  regions  and  will  teminate  in  the  flow 
where  local  Mach  number  is  equal  to  .1.  Then  in  the  region  extending 
to  infinity  behind  the  shock  waves  and  bounded  by  the  stream  lines 
passing  through  the  points  where  the  shock  terminates,  the  flow  is 
rotational  and  with  variable  entropy.  For  convenience,  this  region 
will  be  denoted  by  D2,  and  the  region  outside  D2,  by  Then 

in  D]_  the  following  equations  must  be  satisfied: 

(1)  The  two  equations  of  motion 


UU^  +  VUy  =  -Px/p  I 
UVx  +  Wy  =  -Py/p  j 


(1) 


where  u  and  v  are  the  velocity  components  in  the  x-  and  y-directions, 
p  the  density,  and  p  the  pressure,  and  the  subscripts  x  and  y 
indicate  partial  differentiation  in  the  corresponding  direction. 


(2)  The  equation  of  continuity 


(pu)x  +  (pv)y  =  0 


(2) 

w 


(3)  The  equation  of  state 


P  =  (3) 

where  C]_  is  a  constant  and  7  is  the  ratio  of  specific  heats.  If  the 
subscript  s  is  used  to  indicate  the  stagnation  conditions,  =  Ps^Pg^. 

(4)  The  irrotationality  condition 


Uy  •"  V^  =  0 


(^) 


With  aid  of  equation  (4),  equation  (l)  can  he  integrated  to  give 

(5) 


where  q^  =  and  C2  is  another  constant. 

The  boundary  conditions  are  that  at  infinity  the  velocity  components 
u  and  V  are  equal  to  the  given  values,  and  on  the  solid  boundary  of 
the  airfoil  the  normal  component  of  the  velocity  vanishes. 

In  region  D2  the  flow  becomfes  now  rotational  and  the  entropy  S 
is  no  longer  a  constant.  Here  the  governing  differential  equations  are 
as  follows:  Equation  of  motion  (l)  and  the  equation  of  continuity  (2) 
are  still  valid.  The  equation  of  state  has  to  be  changed  to  the 
following  form: 


__p _ +  =  c 

(7  -  Dp  2  ^^2 


p  = 


p7 


(6) 


where  Cy  is  the  specific  heat  at  constant  volume  and  where  C3  is  a 
constant.  The  energy  equation  is 


uS^  +  vSy  =  0 


(7) 


For  flows  derived  from  isentropic  irrotational  flows  by  a  shock,  there 
is  the  so-called  condition  of  isoenergetlc  flow 


7V  +  al  -  r 
(7  -  1)P  2  ^^2 


(5a) 


where  C2  is  a  constant  throughout  the  flow,  before  and  after  the  shocks. 
This  equation  is  then  identical  to  equation  (5);  it  implies  equation  (7) 
when  equations  (l)  and  (2)  are  satisfied. 

The  boundary  condition  at  the  solid  body  is  the  same  as  before;  that 
is,  the  normal  component  of  velocity  must  be  zero.  The  condition  at 
infinity,  however,  is  no  longer  the  same,  because  if  the  velocity  were 
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constant  at  infinity  the  flow  would  again  be  irrotational.  The  correct 
boundary  condition  in  this  case  is  that  the  pressure  must  be  a  constant 
at  infinity. 

Since  the  flow  in  a  part  of  the  domain  is  rotational^  the  velocity 
potential  cp  obviously  cannot  be-  used  and  it  is  convenient  to  introduce 
in  such  cases  the  stream  function  ^  defined  by 


(8) 


Equation  (7)  indicates  that  the  entropy  S  is  a  constant  along 
each  streamline.  The  entropy  S  can  therefore  be  written  in  the 
following  form: 


S  =  cvf(t) 


(9) 


and  the  equation  of  state  (6)  as 


p  = 


(10) 

» 


With  some  calculation,  the  rotation  of  the  flow  may  be  written  in  terms 
of  as  follows: 


=  Vx  -  Uy 

=  C3P^e^(^h' (t) /(7  -  1) 


(11) 


VARIATIONAL  PRINCIPLES 


Instead  of  studying  the  boundary-value  problems  as  formulated  in 
the  preceding  paragraphs,  it  is  sometimes  more  convenient  to  study  the  » 

associated  variational  problems,  especially  when  the  exact  solutions  of 
the  differential  equations  are  very  difficult  to  obtain.  In  such  cases, 
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approximate  methods  of  solution  such  as  the  Rayleigh-Ritz  method  and 
the  Galerkin  method  will  be  extremely  useful. 

For  irrotational  flows^  two  variational  principles  were  given  by 
Bateman  (reference  8),  one  in  terms  of  the  velocity  potential  and  the 
other  in  tenns  of  the  stream  fimction.  The  first  one  is  more  suitable 
for  stuciying  the  flow  passing  arbitrary  bodies  and  has  been  used  by  the 
senior  author  and  his  associates  in  computing  such  flows  in  refer¬ 
ences  1  to  7.  In  the  present  case,  the  second  principle  which  is 
expressed  in  terms  of  f  should  be  used. 

If  the  pressure  is 


p  =  F(p)  -  pF'(p) 


(12) 


where  F(p)  is  some  function  of  p  and  the  prime  indicates  differentia- 
tion,  the  variational  integral  is 


F(p)  + 


dx  dy 


(13) 


in  which  p  and  ■>|f  are  to  be  varied  independently. 

Equation  (12)  is  a  differential  equation  of  the  Clairaut  type. 
If  p  =  c^^pT' ,  the  solution  of  equation  (12)  is  F(p)  =  cgp  - 
The  integral  Jp  then  becomes 


Cgp 


7-1 


2p 


dx  dy 


(11^) 


The  condition  that  6Jp  =  0  then  leads  to 
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where  the  last  integral  is  taken  on  the  boundary  of  the  domain.  If  this 

line  integral  is  zero,  the  condition  =  0  leads  to  the  following  » 

equations 


7  P 

+  — i - -  c, 

7  -  1  P 


ft/p)x  (♦y/o)y  ■  0 


as  the  Euler  equations  in  the  calculus  of  variation.  The  first  equation 
is  the  Bernoulli  equation  and  the  second  is  the  condition  of  irrotation- 
ality,.  The  condition  of  continuity  is  satisfied  automatically  by  intro¬ 
ducing  the  stream  function  ‘'1^.  Thus  if  the  line  integral  is  equal  to 
zero  the  vanishing  of  the  first  variation  of  J  leads  to  the  desired 
differential  equations.  If  the  line  integral  is  not  zero,  as  in  similar 
cases  discussed  in  references  1,  2,  6,  and  10,  the  value  of  the  integral 
must  be  subtracted  from  5J]_  and  should  be  modified  accordingly. 

It  is  interesting  to  note  that  although  the  variational  integral  (13) 
is  formulated  by  Bateman  for  irrotational  flow  only,  it  is  also  valid 
for  rotational  flow  if  in  equation  (12)  F(p)  is  replaced  by  F(p,S) 
and  is  solved  from  the  correct  equation  of  state,  equation  (6)  or  (lO). 

In  this  case 

r 

F  _  p  ^  =  C3ef(^)p^ 

SP 


•The  solution  of  this  differential  equation  may  be  written  as 

Cqe^^^)p^ 

F(p,Sj  =  Cpp - ^ - 

7-1 


(15) 


and  a  variational  integral  for  the  rotational  flow  after  shock  waves  may 
be  written  as  follows 


(7  -  1)_ 


dx  dy 


(16) 
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This  is  the  same  integral  recently  obtained  by  Lin  and  Rubinov  (refer¬ 
ence  9)-  As  before^  the  condition  6J2  =  0'  leads  to 


i  ds  =  0 

P  on 


The  last  integral  is  again  taken  on  the  boundary  of  the  domain.  If  this 
integral  is  zero,  there  result' 


^  7P 

2p2  (r  -  i)p 


"^3 


7  -  1 


f'(^) 


(17) 


as  the  Euler’s  equations.  The  first  equation  is  the  energy  equation  (5a) 
and  the  second  equation  is  the  rotation  equation  (ll).  With  the  aid  of 
the  second  equation,  the  equations  of  motion  (l)  can  be  derived  from  the 
first  equation  as  follows:  Writing  the  first  equation  in  a  slightly 
different  form. 


+  “{u^  +  v^)  =  C2 

Differentiating  it  with  respect  to  x, 

rc2e^‘(^)p7-2p^  +  -■  7  ^  + 

uu^  +  =  0 


(18) 


4 


10 


NACA  TN  2539 


Since  p  = 


Px/P  =  lC3ef'('l^)p’'-2p^  +  C3e^^'>^^p’'"^f'(t)t, 


Substitution  in  equation  (l8)  results  in 


C3p’'e^(^)f•('^)-  il^x  + 


uu„  +  =  0 

X 


(19) 


As  C3P^e^(^)f' (\lf)/(7  -  1)  =  Vx  -  Uy  and  J  i^x  =  equation  (19)  is 

Wentical  to  the  equation  of  motion  in  the  x-direction.  Similarly^ 
differentiation  of  equation  (17)  with  respect  to  y  results  in  the  other 
equation  of  motion. 


SHOCK  CONDITIONS  AS  A  RESULT  OF  VARIATIONAL  PRINCIPLE 


In  discussing  the  mathematical  formulation  of  the  problem,  the  flow 
as  shown  in  figure  1  may  be  put  into  a  simplified  form  as  shown  in  fig¬ 
ure  2,  where  and  Bp  denote  the  regions  when  the  flow  is  irrota- 

tional  and  rotational,  respectively,  and  C2^_p  is  the  shock  wave  plus 
the  common  boundary.  The  variational  principle  is  that 

6J  =  +  SJg  =  0 


leads  to  the  desired  differential  equations 
Here 


in  the  corresponding  domain. 


Jl 


^1^  ,  frlxg  +  ilyg) 
7-1  2pi 


dx  dy 


(7  -  1)  2P2 


dx  dy  (20) 
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where  as  shown  in  the  above  expressions  is  integrated  over  D]_ 

and  J2  over  D2;  and  P2;’‘t'’2  indicate  that  p  and  i  are 

different  functions  in  the  domains  D]_  and  D2;  and  the  func¬ 
tions  p^^  and  are  to  be  varied  independently.  The 

condition  6J  =  0  therefore  leads  to 


(7  - 1) 


2Pl 


5p 


dx  dy  - 


(^ix/Pijx  (^ly/'^ijy 

dx  dy  + 

ff 

7  7C3ef(^2)p27-l 

(^2x2  +  '^2/) 

IJ  D2 

1 

1 — 1 

1 

2p^  _ 

&P2  dx  dy 


3P/ef('l'2)ffe) 


(7  -  1) 

_  _ L  Bil; 

1  Pi  Sn 


{^2xjP2)x  (^2y/P2)y 


6’^2  ■’* 


r  1  . 

ds  + 


r  1  ^^2  r  /l  ^1 

C2  P2  5n  '^^1-2 Ipi 


St 


6ti  -  ^  6tp  ds  =  0 


Consider  the  line  integrals  first.  On  a  part  of  and  C2  where 

the  boimdary  is  the  surface  of  the  body,  t]^  and  t2  are  equal  to 
the. chosen  constants.  Then  6t]_  =  5t2  =  0,  and  the  first  two  line 

integrals  are  zero.  These  integrals  may  not  be  zero  when  taken  over 
the  other  parts  of  the  boundaries.  In  such  cases,  the  values  of  the 
integrals  must  be  subtracted  from  6J  and  J  must  be  modified  accord¬ 
ingly.  On  the  shock  waves,  if 


(l/pl)(a^|fl/^n)  =  (l/p2)(^'r2/Sn) 


(21) 
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and 


(22) 


the  third  integral  vanishes.  Condition  (21)  requires  the  velocity 
components  tangent  to  a  shock  wave  to  be  the  same  in  passing  through 
the  shock.  Condition  (22)  is  the  continuity  equation  in  crossing  a 
shock  wave.  On  the  common  boundary,  since  it  is  a  streamline, 

=  0  and  the  line  integral  vanishes. 

When  the  line  integrals  are  zero  or  modified  to  be  zero,  the 
following  Euler  equations  are  obtained  from  the  well-known  rules  of  the 
Calculus  of  Variation.  They  are 


+  (^Ix^  +  ^Iv^)  ^ 


(r  -  i)p;l 


2p 


C2  in 


(7  -  1)P2 


2p, 


C2  in  Dg 


(23) 


(2k) 

(25) 


in  D2  (26) 


Equations  (23)  and  (24)  indicate  that  the  energy  in  the  fluid  remains 
unchanged  in  crossing  a  shock  wave.  With  equations  (23)  and  (25)  and 
equations  (24)  and  (26),  the  equation  of  motion  can  be  derived;  that  is, 
the  momentum  of  the  flow  also  remains  the  same  in  crossing  a  shock.  The 
continuity  equation  is  already  given  as  equation  (22).  Therefore,  from 
the  condition  of  6J  =  0,  all  the  conditions  in  crossing  a  shock  wave 
are  obtained. 
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SUBSONIC  ANT)  TRANSONIC  FLOWS  PAST  A  CIRCULAR  CYLINDER. 
BEFORE  OCCURRENCE  OF  SHOCK  WAVES 


As  a  first  step  to  the  solution  of  transonic  flows  with  shock  waves ^ 
the  flow  before  the  occurrence  of  shock  waves  will  be  considered.  The 
variational  integral  J  contains  two  variables^  p  and  Numerical 

calculation  may  be  carried  out  by  following  the  Rayleigh-Ritz  procedure. 
Unlike  the  cases  considered  in  references  1  to  6  where  the  velocity 
potential  cp  is  used,  the  labor  involved  in  the  computation  by  taking 
and  p  as  variables  in  the  Raylei^-Ritz  method  becomes  excessively 
large.  Instead,  Galerkin’s  method  was  found  to  be  much  more  simple 
in  this  case.  The  application  of  Galerkin^s  method  to  compressible- 
flow  problems  has  been  discussed  in  detail  in  references  3  and  4  and 
can  be  briefly  outlined  as  follows:  First  the  variables  ijr  and  p 
are  written  in  the  form  of  series  which  satisfy  the  boundary  conditions 
but  with  undetermined  parameters,  such  as 


^  +  51  Z  V 

m  n 


> 


m  n 


(27) 


where  A^  and  B^^  are  the  undetermined  parameters.  Then  A^^ 
and  Bjjin  are  determined  from  the  conditions 


With  and  B^^^  determined,  equations  (27)  give  an  approximate  solu¬ 

tion  to  the  problem. 


It  was  found,  however,  that  near  the  limiting  Mach  number  where 
shock  waves  are  about  to  occur,  the  p-series  converges  rather  slowly  and 
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many  parameters  have  to  be  used.  The  amount  of  work  involved  in  the 
integration  and  in  the  solution  of  the  resulting  simultaneous  equations 
with  many  parameters  again  becomes  very  great.  A  modified  Galerkin 
method  has  been  tried  and  it  was  found  possible  to  reduce  much  of  the 
numerical  work.  The  method  is  essentially  as  follows:  Instead  of 
determining  by  Galerkin ^s  method,  these  parameters  can  be  solved 

in  terms  of  A^n  from  the  Bernoulli  equation: 


("2^^  "  7^  -  2  ^  v)  =  ° 

by  the  method  of  equal  coefficients.  The  parameters  then  become 

the  only  unknowns  and  they  are  determined  by  the  Galerkin  method. 

In  the  case  of  two-dimensional  flow  past  a  circular  cylinder'  with 
unit  radius  the  boundary  conditions  are  at  r  =  <»; 


P  =  Po 

fyelPo^)^  +  (^/po)^  =  u2 


and  at  r  =  1 


^0  =  0 


Consider  the  case  where  the  circulation  is  zero.  The  flow  must  then 
be  symmetrical  with  respect  to  both  the  x-  and  y-axis.  The  boundary 
conditions  and  the  conditions  of  symmetry  are  satisfied  if  ^  and  p 
are  assumed  in  series  of  the  following  forms: 


-.SI;  *■■(*- 


sin  n9 


(31) 


and 


P 


00  00 


m=l  n=0 


cos  2n0 
r^m 


(32) 
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Taking  two  parameters  in  ilf,  equation  (31)  "becomes 

—  -t-, 

(r  -  i)  sin  0  +  Aii(i  -  sin  0  +  Ai3^1  -  sin  30 


t  =  P^U 


Substituting  equation  (33)  into  equation  (30)  and  equating  the  coeffi 
cients  of  ^  and  cos  n9^  there  result 


'10 


®11  = 


B. 


In 


B. 


B. 


= 


B. 


0 

L 

4 

%  -  IV 

2Mo 

i 

1 

1— 1 

1 

13 

0  if  n  ^ 

1 

M  2  -  1 

2 

0 

1 

1 — 1 

1 

OJ 

0 

S 

2 

1 

Mo^ 

V  -  1! 

2 

1 

V  -  1 

2 

'  2  2 

All  ^l^SAij-' 


l-kAii  *  *AiiAi3  -  ^13^  .  BiiBi^H 


\  1^  2 

13  ■  ^"iri3  j  ■  2 


■SA,.  +  2A_A. 


Bii-H 


(‘*''13")  - 


(33) 
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B 


1  1  -a  2tT 

2^  "  Mo2  -  1  ?  12  H 


=0  if  n  i  5 

2n 


B 


30  ■  M  2  .  1 


2M, 


(^11  ■  -^11^  "  3Ai3  ^  ^ 


B 


31  Mo2  -  1 


^o^(--^ll  +  ^11^  +  3Aj^2  "  ^11^3)  "  (2^12^20  ■*■  ®11®22 


®12®21  ■*■  ^12^23)  ^  6  ^  i  ®11®12^) 


B 


32  -  1 


'  '  (®lAl  "  ®11®23  "  "®12®20  *  »12®21.)‘'  * 


-B  ^'l 
2  12  / 


B 


33  M/  -  1 


-3M/A13B  -  .  (>• 


)H  .  ^(. 


V" 


3BuBi2^) 


B 


3^  Mo2  _  1 


-( 


«11®23  *  *12*22)“  *  *ll"*12 


B 


35  -,^2.1 


-  B  B  .  +  B  B 

\  11  2k  12  23 


W  .  i2L^  B  B  * 


B 


36  "  M^2  -  1 


-*12*24*  *  ^3-^  B^^s 


3 
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B.  =  0  if  n  ^  7 


1^0  m2 


^  +  »li^)  -  "(=20^  t  i^21^  *  +  Ib^^; 


2^24^  +  ■*■  B12B32)  ■*■ 


*11^^20  1*11^^22  * 


5211*12*21  +  5®11®12*23  *  ^2^0  +  ^12%k)  *  ^Hl^d*!!**  + 


3-0  2-r,  2  3-0  V 

2^11  ®12  +  H®12 


1  (2  \  f 

^  2  _  2.  ^  ^^^ll'^is)  ■  ^(^®20®21  ®21®22  ''■  ®22®23 


20-°21  21  22  22  23 


®23®2l|  ^^11^30  ®11®32  ^12^31  ■*■  ^12^33)  '2  (^®12^®21 

f%l^®23  ■*■  ®11®12®20  ®11®12®23  2^11^12^24  2®12^®21 


t®12^®23)  ^ 


/2B^  3b  +  Ib  B  3 
\  11  12  2  11  12 


’^2  “  1  "3^0%iAi3  -  +  ^\o\2  ^21^3  +  ®22®24  ^  ®11®31 

%1®33  ^  2®12®30  \2®34)  ^  ^11%2  i\i\k 

®11®12®21  +  |®11®12®23  4^12%2)  ■(?'  + 


®11®12®21  +  |®11®12®23  4^12%2)  '"" 


-B, + 


”R  ^ 

®11  \2 
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®'»3  "  *  *21=22  "  *21*24  ♦  **11*32  "  **11*34  " 

**12*31 "  **12*35)  *  ^^4^(Ki%i  ♦  1*11"*; 


'23 


*11*12*20  '*'  2*11*12*22  '*'  2*11*12*24  ■*■  5'*12  *21  2*12  *23) 


+  +  B^  ^B^^  +  ^11^35  ^12^33 


*'‘*t  ”  m^2''.  iE(i*22^  *  ^*20*24  ■*  ^21*23  * 

*12^36)  *  ^*-H^{i*llS2  *  i®ll\4  +  1*11*12*21  *  1*11*12*23 

4*12**20  *  4*12**24)  -  -  *)  feBii-^  ♦  +  ^J) 


B 


1^5 


^  ^22^23  ^  \i"34  "ii"36  ^  ^2^33) 

2'  — ■(^11^B23  +  |Bj^^B^2®22  2®ll®12®2l).  +  li:®12^®2l) 


(7  1)|;.(7  ■--21(1b,,33^,  .  1b,^b,23 


46  m-2 


Mo^  -  1 


-*(4*23"  "  *22*24  *  *11*35  *  *12*34)  *  H^(4*1i"*24  * 


i®ll®12®23  F 


i»i2%3) 


(7  +  1)7(7  -  2)/3t,  2t,  2\ 

\f\i  ®12  j 


W 


B 


4? 


3  3^  +  (7  +  IJi/Ie  B  B  ,  +  ^B  ^B  \  + 

12^35;  2  \2  11  12  24  4  12  23/ 
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B 


48 


Mo^  -  1 


-“(Ki."  "  Vie)  *  * 


) 


(7  +  1)7(7  -  2) 

25 


(34) 


=  0  if  n  ^  9 


where  H  =  ^9 - 7  (/ ..,±  j-)  ~  For.  given  free-stream  Mach  nimher  and  7, 

2  2(7  -  1) 

the  Ei^’s  are  now  expressed  in  terms  of  and 

numerical  example^  12  parameters  are  taken  in  the  p-series  which  seems 
to  give  satisfactory  convergence  at  hi^  Mach  numbers. 

The  condition  of  irrotationality  is  now  the  only  differential 
equation  remaining  to  he  satisfied.  In  terms  of  if  and  polar  coordinates^ 
equation  (4)  becomes 

Pr^p  +  -  P(Ailf)  =  0  (35) 


4 


where 


At  = 


"rr 


tee 

i2 


The  Galerkin  equation  for  determining  the  parameters  is: 


r=oo  n  e=2)t 


r=l'^  9=0 


Prt] 


+  Pete 

“3“ 


-  p(At) 


r) 


8t 


r  de  dr  =  0 


(36) 
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With  the  two  parameters  A^ 


and  A-|_2 


in 


equation  (36)  becomes: 


0  d0  dr  =  0 


(37) 


30  d0  dr  =  0 


(38) 


Carrying  out  the  integration  and  taking  7  =  1.405  the  following 
equations  are  obtained  for  different  values  of  Mq.  The  parameters 

taken  in  the  p-series  are  ^20f  ^21>  ®22>  ®23^  ®30^ 

®31^  ®32^  \l^  \2‘  these  there  result 

at  Mq  =  0.3: 

0.  7914i66A3_]^  +  0.0i48886A3_3  =  -0.0529632  +  O.OOll372Aj^j^A^2  " 

0.0290381Aj^3^2  _  o.0353153A^3^  + 

0.  028i407Aii^Aj^2  -  0.  0368037A^^A^2^  - 
0.026577Aii3  -  0.  0537032A3^33 

-0.04i6032A^3_  +  2.  1113662A;i^3  =  0.0565330  -  0.0524788A^^A^^  - 

0.0049465A^^2  _  o.0516919Aj^2^  - 

0.0k26662A^j_^A^^  +  0.  0614285A^^A^2^  + 

0.0182144a^j^3  _  o.0776527A^23 
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At  Mq  =  0.4: 

0.  8939718Aj_3^  -  0.0293045A;L3  =  -0.1386222  -  0. 0363370A^^A^^  - 

0.  0208205Aj^3_2  _  0. 159415  1A;L3^  + 

0. 0557317Aii^A^2  "  0. 1400827A^^A^22 
0.  0837686A3^3^^  -  0.  i685111A^2^ 

-0.  II65173A11  +  2.21684o8A3_3  =  0.1313933  -  0.  i4623miiA;L3  - 

0.1092484Aj^j^2  -  0.20k6'je9k-^-^  - 
0.053295Aii2a^3  +  0. 3247757AilA3^3^ 

0.059l8llA3_j^3  _  0.  1892021A3^2^ 

At  Mq  =  0.48: 

1.  0588773A^3^  -  0.0792806Aj^^  =  -0.3375697  -  0.0768956a^^A^2  - 

0.  030879A3^3^^  -  0.7005271A^2^  + 

0.  0580232A3^3^2^j^2  +  0.  362238A^^A^^2  _ 

0.  0786318A3_^3  _  0.  1904891A^2^ 

-0.2295625A.j,  ]_  +  2.215112A3^3  =  0.2795736  -  0.7384032Aj^3^A^2  - 

0.227209IA112  _  0. 9122836A132  + 
0.210928Aii2a^3  +  1.  067842A;j^3^A3_32  + 

0.  1269145A;li^  -  0.  i446766Aj_33 
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At  Mq  =  0.5: 

1.  058801<-5Aii  -  0.  0943961)-A]^3  =  -0.l^35^569  -  0. 074l4-038AiiAi3  - 

0. 0687631Aii^  -  1.0282597Ai3^  - 
0.  0279679Aii^A3_3  +  0,  390i^425Aj_^Aj^2^  + 

1.  236569A3_j^3  _  0.579232aAi33 

-0.2586206A11  +  2. 15l^37l^9Ai3  =  0.3^^71313  -  0.9668059A^3^A^2  - 

0.289680lt.Aj^l2  _  i.2930194Ai3^  + 

0.4529877A3_;L^A^2  +  1. 7020295A^]_A^2^  + 

0.136i617Aii3  -  0.3869221^3^2'^ 


These  equations  can  he  solved  by  the  method  of  successive  approxi¬ 
mation  as  outlined  in  references  1  and  2.  The  parameters  at  various 
Mach  numbers  and  the  computed  maximum  velocities  at  these  Mach  numbers 
are  shown  in  table  I.  These  values  show  good  agreement  with  those 
obtained  by  other  approximate  methods.  At  Mq  =  0.5  the  method  of 
successive  approximation  fails ^  and  it  was  found  by  graphical  method 
that  these  equations  do  not  have  a  solution.  This  agrees  with  a  previous 
investigation  as  reported  in  reference  7  where  the  velocity  potential  cp 
was  used.  This  indicates  that  the  flow  is  no  longer  irrotational  and 
shock  waves  have  probably  appeared.  The  limiting  free-stream  Mach 
number  is  between  0.48  and  0.5«  This  value  is  probably  too  high  because 
only  two  parameters  in  are  used.  This  fact  has  been  discussed  in 
reference  7* 

To  get  an  idea  of  the  convergence  of  the  and  p-series,  the  values 
of  maximum  velocities  and  p  at  Mq  =  0.48  are  computed  and  listed  in 

tables  II,  III,  and  IV.  The  p-series  does  not  appear  to  be  convergent 
when  the  parameters  are  added  one  by  one,  but  if  these  parameters 

are  added  by  the  group  the  series  becomes  convergent.  This  can  be  seen 
by  comparing  the  values  of  p  and  the  maximum  velocities  with  various 
numbers  of  as  given  in  tables  II  and  III. 
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TRAWSOSriC  FLOW  PAST  A  CIRCULAR  CYLINDER  WITH  SHOCK  WAVES 


When  the  free-stream  Mach  numher  reaches  0. 5  there  is  no  solution 
to  the  equations  obtained  by  assuming  irrotational  flow  in  the  entire 
domain.  Shock  waves  probably  have  occiarred.  Thus  the  flow  is  irrota¬ 
tional  only  in  the  region  before  the  shock  waves  or  D]_  and  becomes 

rotational  in  the  region  behind  the  shock  waves  D2.  The  approximate 

solution  according  to  the  modified  Galerkin  method  can  be  carried  out 
as  follows: 

(1)  The  locations  and  the  obliquity  of  the  shock  waves  are  first 
assumed  from  the  best  evidences  available. 

(2)  When  the  assumed  shockwaves  are  Inserted  in  the  flow,  the 
domains  Dq  and  D2  are  thus  defined.  In  the  region  Dq,  the  flow 

is  still  irrotational.  The  boundary  conditions  are  that  the  velocity 
at  infinity  is  equal  to  the  undisturbed  velocity  and  on  the  surface  of 
the  body  must  be  equal  to  a  chosen  constant.  The  values  of  i|fq 

and  pq  may  be  assumed  in  the  same  form  as  in  the  case  before  the  shock 
waves  occur.  The  parameters  A^j^  and  are  determined  from  the 

modified  Galerkin  method.  The  only  difference  between  the  present  case 
and  the  previous  case  is  that  in  this  case  the  Integrals  are  to  be 
extended  only  in  the  region  Dq  Instead  of  the  entire  flow  region. 

(3)  The  entropy  distribution  S  =  c^f(\|()  in  the  flow  after  the 
shock  is  to  be  assumed  next.  From  the  flow  conditions  before  the  occur¬ 
rence  of  shock  waves,  a  good  approximation  of  fCi)  may  be  found  by 
rough  preliminary  calculations. 

(h)  The  "boundary  condition  at  infinity  of  the  flow  after  the  shock 
is  that  p  is  equal  to  the  constant  pressure  p^.  With  f(t)  assumed^ 

computed  from  equation  (lO)  and  the  velocity  distribution 
from  equation  (5a). 

(5)  With  the  above  boundary  conditions  and  the  condition  that  =  ^2 

on  the  assumed  shock  waves  and  common  streamlines’*  of  and  D2^  ^2 

and  p2  are  then  assumed  in  the  forms  of  series  with  undetermined  param¬ 
eters.  These  parameters  can  again  be  determined  by  the  modified  Galerkin 
method. 

(6)  After  these  undetermined  constants  in  and  ^2  Q-re  computed^ 

the -values  of  (l/pi)  and  (1/P2)  (aitfg/Sn)  at  the  shock  waves 
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can  "be  foimd  and  is  calculated.  If  is  equal 

to  ^l/P2j  (^^2/^^)  foimd  is  not  the  same  as  f(l^)  assumed, 

as  is  generally  the  case,  adjust  the  locations  and  the  obliquity  of  the 
shock  waves.  Then  based  on  the  new  position  and  shape  of  the  shock  waves, 
assume  a  new  f(^)  and  repeat  the  computations  agsin  until  these  values 
check  with  each  other  in  two  consecutive  cycles. 

A  numerical  example  will  now  be  carried  out.  With  and  p]_  as 

assumed  in  equations  (3I)  and  (32)  and  by  taking  2  parameters  An 
and  A]_3  in  equation  (31)  and  12  parameters  in  equation  (32),  it  was 
found  that  there  does  not  exist  a  solution  at  Mq  =  0.5*  In  refer¬ 
ence  11  it  is  shown  that  when  the  shock  waves  first  occur  they  are 
almost  normal  shocks  located  near  the  point  of  maximum  velocity.  As 
a  first  approximation,  assimie  that  at  Mq  =  0. 5  there  are  two  normal 
shock  waves  located  at  9  =  n/2  and  3^/^  from  r  =  1  to  r  =  1.5. 

The  separating  streamlines  between  the  regions  D]_  and  D2  are  assumed 
to  be  the  same  as  the  streamlines  dn  the  flow  at  Mq  =  0.48  passing 
through  the  points  r  =  1.5  and  9  =  jt/2,  3^2* 

To  determine  the  parameters  and  integrals  (37)  and  (38) 

are  to  be  extended  in  the  above  assumed  region.  This  can  be  done  by 
taking  first  the  integration  over  the  entire  flow  region  and  then  sub¬ 
tracting  from  the  resulting  values  the  values  of  these  integrals  evalu¬ 
ated  in  the  region  D2.  The  region  D2  consists  of  two  areas,  namely, 

a  semicircular  ring  with  r  from  1  to  1.5  and  9  from  it/2  to  -7t/2  and 
a  symmetrical  tail-shaped  area.  The  dominating  parts  of  the  values  of 
the  integrals  happen  to  be  in  the  semicircular-ring  area  where  the 
integration  can  be  carried  out  analytically.  In  the  tail-shaped  area, 
analytical  integration  becomes  rather  laborious  and  numerical  integration 
has  been  used.  After  subtraction,  these  two  equations  are 

2.0587561^Al3_  -  0.2021789A13  =  -0.8773276  -  O.1912757A11A13  - 

0.  1931520A3^3^2  _  2.01076aA3^32  + 

0. 3*^9061f8A^3^2^j_3  +  0.89i^8537A3^^A^^2  _ 

O.OO87952A11  +  9. 2190418A133 

J| 
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-0.335855i+Aii  +  5.2963361A13  =  0.76527^8  -  1. 7476647AiiAi3  - 

0.1i.62l^381^A]_3_2  _  2.088238Aj^3^  + 
0.7361791Aii2a^3  +  2. 9033634a^^A^32  + 
0.3571j436A3^3^3  _  ,x919k2hk-^^^ 

While  a  solution  does  not  exist  "before  the  integrals  in  D2  sub¬ 

tracted,  solution  of  these  equations  now  does  exist.  Solving,  there 
result: 

All  -0-^37  A3_3  =  0.100 

With  these  values  of  A-^i  and  Aj^3,  the  point  at  B  =  n/2  with  a 

unit  local  Mach  number  is  computed  to  be  at  r  =  1,6.  This  indicates 
that  the  shock  waves  probably  would  terminate  at;  a  value  of  r  between 
1.5  and  1.6,  and  the  region  D]_  as  assumed  is  a  good  approximation. 

The  boundary  conditions  for  '^2  P2  follows: 

At  r  =  60,  p  =  Po  and  therefore 

p/  .  Sa  (39) 

At  r  =  since  is  taken  as  zero^ 

t2  =  0  (40) 

To  insure  the  vanishing  of  normal  velocity, 

=  0  (4l) 


On  the  shock  waves  and  the  common  boundary, 

a 


^2  =  '^1 


(42) 
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Conditions  (39)^  (^O),  and  (4l)  are  satisfied  if  '^2  and  Pg  are 
assumed  in  the  following  forms: 


and 


P2 


1/7 


00 


00 


*ZZb' 

r=oo  m=l  n=0 


mn 


cos  2n0 
^2m 


m 


Note  that  becomes  a  function  of  y  =  r  sin  0  at  r  =  00.  The 

constants  are  to  be  determined  by  the  collation  method  in  terms 

of  that  at  a  number  of  points  on  the  shock  waves 

and  the  common  boundary.  With  ^{^2)  assumed,  may  be  solved 

in  terms  of  A’^^  from  equation  (24)  and  the  parameters  may 

be  calculated  from  equation  (26)  by  using  Galerkin's  method.  A  numeri¬ 
cal  example  has  been  tried.  It  was  found,  however,  that  unless  some 
modern  hi^-speed  computing  machine  is  used  the  computation  becomes 
very  lengthjr. 


New  York  University 

New  York,  N.  Y. ,  March  1,  I95I 
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TABLE  I 

PARAMETERS  AIO)  COMPUTED  MAXIMUM  VELOCITIES 
AT  VARIOUS  MACH  NUMBERS 


M. 


0.3 


0.4 


0.47 


0.48 


11 


-0. 0677 


-0.1537 


-0.296 


0.316 


■13 


0. 0255 


0. o4o4 


0.088 


0.087 


at 


r  =  1 


e  = 


at 


r  =  1 


e  =  90^ 


0. 8254 


0.  6972 


0. 3705 


2. 1972 


2.2859 


3.4870 


0.2680 


4.5525 
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TABLE  II 

CONVEEGMCE  OF  p-SERIES  AT  =  0.48 


Parameters 
taken  in 
p  series 

Bll 

B11B12 

B11B12 

®20 

®11®12 

B20B21 

B11B12 

B20B21 

B22B23 

B11B12 

®20®21 

B22B23 

B30B31 

B32 

®11®12 

B20B2I 

B22B23 

B30B31 

B32B40 

B4iB1^2 

^11 

-0. 0988 

-0.0836 

-0.490 

-0.390 

-0. 3950 

-0. 3547 

-0.316 

■^13 

0,0777 

0,0784 

0.120 

0.091' 

0.0960 

0,1110 

0,074 

r  =  1 

^  at 

0  =  90° 

0. 6943 

0. 6522 

0. 14l6 

0.2931 

0.2686 

0.2788 

^2680 

q  ■  ^ 

3  at 

^  0  =  90° 

2.3765 

2.5698 

5.5104 

3.5416 

3.7902 

4.5503 

4.5525 

a  |>n 
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TABLE  IV 

♦ 

CONVERGENCE  OF  t  SERIES  WITH 
12  PARAMETERS  IN  p -SERIES  AT  Mq  =  OA8 


Parameters  taken 
in  't-series 


P  ^  ^ 

^  at 

^  0  =  90^ 


r  =  1 


e  =  90^ 


1 — 1 
t — 1 

c 

^11^13 

0. 1007 

0.2680 

13.519i^ 

4.5525 

NACA^ 
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breakdown  of  the  irrotational  flow.  By  allowing  breakdown  of  the  irrotatlonal  flow.  By  allowing 

shock  waves  to  occur,  solution  again  exists.  shock  waves  to  occur,  solution  again  ffldsts. 


